//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <stdint.h>

#include <vector>
#include <string>
#include <map>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#ifdef QUAD_VALUE
  typedef long double value_t;

  inline value_t vpow(const value_t& x, const value_t& y)
  {
    return powl (x, y);
  }

  inline value_t vpowi(const value_t& x, int y)
  {
    return powl (x, y);
  }

  inline value_t vsqrt(const value_t& x)
  {
    return sqrtl (x);
  }

  inline long v2si(const value_t& x)
  {
    return (long)x;
  }

  inline double v2d(const value_t& x)
  {
    return (double)x;
  }

  inline long double v2dd(const value_t& x)
  {
    return x;
  }

  inline value_t dd2v(const long double& x)
  {
    return x;
  }
#elif defined(GMP_VALUE)
  #include <gmpxx.h>

  typedef mpf_class value_t;

  inline value_t vpow(const value_t& x, const value_t& y)
  {
    return (double)powl (x.get_d(), y.get_d());
  }

  inline value_t vpowi(const value_t& x, int n)
  {
    if ( n == 0 )
      {
        return 1;
      }
    else if ( n < 0 )
      {
        value_t ret = x;
        for ( int i = n; i < -1; i++ )
            ret = ret * x;
        return 1.0 / ret;
      }
    else
      {
        value_t ret = x;
        for ( int i = 1; i < n; i++ )
            ret = ret * x;
        return ret;
      }
  }

  inline value_t vsqrt(const value_t& x)
  {
    return sqrt (x);
  }

  inline long v2si(const value_t& x)
  {
    return x.get_si();
  }

  inline double v2d(const value_t& x)
  {
    return x.get_d();
  }

  inline long double v2dd(const value_t& x)
  {
    return x.get_d();
  }

  inline value_t dd2v(const long double& x)
  {
    double ret = x;
    return ret;
  }
#else
  typedef double value_t;

  inline value_t vpow(const value_t& x, const value_t& y)
  {
    return pow (x, y);
  }

  inline value_t vpowi(const value_t& x, int y)
  {
    return pow (x, y);
  }

  inline value_t vsqrt(const value_t& x)
  {
    return sqrt (x);
  }

  inline long v2si(const value_t& x)
  {
    return (long)x;
  }

  inline double v2d(const value_t& x)
  {
    return x;
  }

  inline long double v2dd(const value_t& x)
  {
    return x;
  }

  inline value_t dd2v(const long double& x)
  {
    return x;
  }
#endif

//#define UNCORRELATED
//#define FREQANAL

#define SQRT2 1.414213562373095048801688724209698

#ifdef UNCORRELATED
  #undef CUMULATION
  #define CUMULATION 1
#else
  #ifndef CUMULATION
    #define CUMULATION 5
  #endif
#endif

#define __PERC_INLINE __attribute__((always_inline))
//#define __PERC_INLINE

extern double EXPONENT;
extern double BASEPOW;
extern int LAYERSTART, LAYEREND;

#ifndef DISTPOWER
  #define DISTPOWER 2.0
#endif

#define CHUNK_LEVEL 5
#define CHUNK_SIZE (1 << CHUNK_LEVEL)

#define ALPHA 1.0
#define BASELINECORRECTION 1
#define RND_VARIANCE (1.0/12.0)

#ifdef UNCORRELATED
  #undef LAYERCOUNT
  #define LAYERCOUNT 0
#else
  #ifndef LAYERCOUNT
    #define LAYERCOUNT 30
  #endif
  #ifdef FREQANAL
    #error "UNCORRELATED cant be on with FREQANAL"
  #endif
#endif

#ifndef FIELD_EXP
  #define FIELD_EXP 12
#endif
#define FIELD_SIZE (1 << FIELD_EXP)

#define RANDOM_DEVICE "/dev/urandom"

typedef long long unsigned int big_t;

#define RNDF_BLOCK_SIZE (7*10000)

//#define RNDF_MMAP

class RNDF
{
public:
  RNDF ();
  ~RNDF ();

  void seed (const gsl_rng_type *type, unsigned long int n, const char *binpath = "rnd.bin");

  inline double rnd () __PERC_INLINE { return next() - 0.5; }
  inline double rnd1 () __PERC_INLINE { return next(); }

private:
  gsl_rng *rng;
  union {
    unsigned char bits[8];
    double dbl;
  } b;

#ifndef RNDF_MMAP
  unsigned char *dptr;
  FILE *fp;
  size_t flen;

  unsigned char data[RNDF_BLOCK_SIZE];
  unsigned char block1[RNDF_BLOCK_SIZE], block2[RNDF_BLOCK_SIZE];
#else
  int memfd;
  size_t flen;
  unsigned char *mem;
  unsigned char data[RNDF_BLOCK_SIZE];
  unsigned char *dptr;
#endif

  void reload();
  double next();
};

extern RNDF rndf;

extern value_t calc_cell();

static inline value_t dist2hf (const value_t& dist2) __PERC_INLINE;
static inline value_t dist2hf (const value_t& dist2)
{
#ifdef HF_FLAT
  if ( dist2 < 1.5 * 1.5 )
      return 1;
  else
      return 0;
#elif defined(HF_CONE)
  if ( dist2 < 1.5 * 1.5 )
      return (1.0 - vsqrt(dist2) / 1.5);
  else
      return 0;
#elif defined(HF_FAST)
  if ( dist2 < 1.5 * 1.5 )
      return 1.5 * 1.5 - dist2;
  else
      return 0;
#else
  #error "invalid or undefined HF"
#endif
}

struct vect2d
{
  double x, y;

  inline vect2d () { }
  inline vect2d (double x, double y) : x(x), y(y) { }

  friend inline vect2d operator+ (const vect2d& l, const vect2d& r) { return vect2d(l.x + r.x, l.y + r.y); }
  friend inline vect2d operator- (const vect2d& l, const vect2d& r) { return vect2d(l.x - r.x, l.y - r.y); }
  friend inline double operator* (const vect2d& l, const vect2d& r) { return l.x * r.x + l.y * r.y; }

  friend inline vect2d operator* (const vect2d& v, double d) { return vect2d(v.x * d, v.y * d); }
  friend inline vect2d operator* (double d, const vect2d& v) { return vect2d(v.x * d, v.y * d); }

  friend inline vect2d operator/ (const vect2d& v, double d) { return vect2d(v.x / d, v.y / d); }

  inline const vect2d& operator+= (const vect2d& v) { x += v.x; y += v.y; return *this; }
  inline const vect2d& operator*= (double d) { x *= d; y *= d; return *this; }
  inline const vect2d& operator-= (const vect2d& v) { x += v.x; y += v.y; return *this; }

  inline double length () const { return sqrt(*this * *this); }
  inline vect2d norm () const { return *this / length(); }
};

class MemPool
{
#ifdef PERC_MEMALIGN
  #define MEMPOOL_LINK_SIZE PERC_MEMALIGN
#else
  #define MEMPOOL_LINK_SIZE sizeof(link_t)
#endif
public:
  MemPool (size_t size, size_t grow);
  ~MemPool ();

  inline void * get () __PERC_INLINE
  {
    if ( begin == end )
      {
        if ( !current->next )
            extend();
        current = current->next;
        begin = ((char *)current) + MEMPOOL_LINK_SIZE;
        end = begin + size * grow;
      }
    void *ret = begin;
    begin += size;
    return ret;
  }

  inline void release (void *mem) __PERC_INLINE
  {
  }

  void reset ();

private:
  struct link_t
  {
    link_t *next;
  };

  size_t size;
  size_t grow;
  link_t *chunks, *last, *current;

  char *begin;
  char *end;

  void extend ();
};


class Chunk
{
public:
  Chunk (int xc, int yc) : xc(xc), yc(yc)
  {
    memset (vals, 0, sizeof (vals));
  }

  double get (int x, int y) __PERC_INLINE
  {
    if ( !vals[x][y].i )
        vals[x][y].d = rndf.rnd();
    return vals[x][y].d;
  }

  template <class RNG>
  inline double get (int x, int y, const RNG& rng)
  {
    if ( !vals[x][y].i )
        vals[x][y].d = rng();
    return vals[x][y].d;
  }

  static inline void * operator new (size_t size)
  {
    return pool->get();
  }

  static inline void operator delete (void *mem, size_t)
  {
    if ( mem )
        pool->release (mem);
  }

  static void reset () { pool->reset(); }

  inline int get_xc () const __PERC_INLINE { return xc; }
  inline int get_yc () const __PERC_INLINE { return yc; }
  Chunk *next;

private:
  union val_t {
    double d;
    uint64_t i;
  };
  val_t vals[CHUNK_SIZE][CHUNK_SIZE];

  int xc, yc;

  static MemPool *pool;
};

class Database
{
private:
  int m_chunkcount;
  Chunk *chunks;
  Chunk **data;

public:
  Database (int m_chunkcount) : m_chunkcount(m_chunkcount)
  {
    chunks = 0;
    data = (Chunk **)calloc (m_chunkcount * m_chunkcount, sizeof (Chunk *));
    reset ();
  }

  inline void reset() __PERC_INLINE
  {
    while ( chunks )
      {
        data[(chunks->get_xc() >> CHUNK_LEVEL) + (chunks->get_yc() >> CHUNK_LEVEL) * m_chunkcount] = 0;
        chunks = chunks->next;
      }
  }

  inline Chunk *get_chunk (int64_t cx, int64_t cy) __PERC_INLINE
  {
    //assert ( cx >= 0 && cy >= 0 && cx < m_chunkcount && cy < m_chunkcount );
    Chunk *c = data[cy * m_chunkcount + cx];
    if ( !c )
      {
        c = new Chunk (cx << CHUNK_LEVEL, cy << CHUNK_LEVEL);
        data[cy * m_chunkcount + cx] = c;
        c->next = chunks;
        chunks = c;
      }
    return c;
  }

  inline double get_by_layer (int lx, int ly) __PERC_INLINE
  {
    int cx = (lx >> CHUNK_LEVEL);
    int cy = (ly >> CHUNK_LEVEL);
    Chunk *c = get_chunk (cx, cy);
    return c->get(lx - (cx << CHUNK_LEVEL), ly - (cy << CHUNK_LEVEL));
  }

  template <class RNG>
  inline double get_by_layer (int lx, int ly, const RNG& rng)
  {
    int cx = (lx >> CHUNK_LEVEL);
    int cy = (ly >> CHUNK_LEVEL);
    Chunk *c = get_chunk (cx, cy);
    return c->get(lx - (cx << CHUNK_LEVEL), ly - (cy << CHUNK_LEVEL), rng);
  }
};

class BASERNG
{
  std::vector<value_t> heights;

public:
  void init (int base0, int base1)
  {
    heights.clear();
    for ( int level = base1; level >= base0; level-- )
      {
        for ( int j = 0; j < CUMULATION; j++ )
          {
            value_t height;
            if ( EXPONENT < 0 )
                height = vsqrt (ALPHA * (vpow (BASEPOW, -EXPONENT) - 1)) * vpow (BASEPOW, -EXPONENT * (level / 2.0));
            else
                height = vsqrt (ALPHA * (1 - vpow (BASEPOW, -EXPONENT))) * vpow (BASEPOW, -EXPONENT * (level / 2.0));
            printf ("BASERNG: level=%d, height=%g, blocksize=%g\n", level, v2d(height), v2d(vpow (BASEPOW, level)));
            heights.push_back (height);
          }
      }
  }

  inline double operator() () const __PERC_INLINE
  {
    value_t ret = 0;
    for ( size_t i = 0; i < heights.size(); i++ )
        ret += calc_cell() * heights[i];
    return v2d(ret);
  }
};

class Layer
{
private:
  int level;
  value_t m_blocksize, m_1dbs;
  int m_chunkcount;
  Database db;
#ifndef UNCORRELATED
  value_t moffsetx, moffsety;
  value_t cosphi, sinphi;
  value_t height;
#endif

public:
  Layer (int level) :
      level(level),
      m_blocksize(vpowi (BASEPOW, level)),
      m_1dbs(1.0 / m_blocksize),
      m_chunkcount(v2si(1 + (5 + SQRT2 * FIELD_SIZE / m_blocksize) / CHUNK_SIZE)),
      db(m_chunkcount)
  {
#ifndef UNCORRELATED
#ifdef FORCE_HEIGHT
    height = 1;
#else
    if ( EXPONENT == 0 )
        height = 1;
    else if ( EXPONENT < 0 )
        height = vsqrt ( ALPHA * (vpow (BASEPOW, -EXPONENT) - 1)) * vpow (BASEPOW, -EXPONENT * (level / 2.0));
    else
        height = vsqrt ( ALPHA * (1 - vpow (BASEPOW, -EXPONENT))) * vpow (BASEPOW, -EXPONENT * (level / 2.0));
#endif
    printf ("Level=%d, blocksize=%g, chunkcount=%d, height=%f\n",
        level, v2d(m_blocksize), m_chunkcount, v2d(height));
#else
    printf ("Level=%d, blocksize=%g, chunkcount=%d\n",
        level, v2d(m_blocksize), m_chunkcount);
#endif
    reset ();
  }

  value_t value (double dx, double dy) __PERC_INLINE
  {
#ifdef UNCORRELATED
    return db.get_by_layer (dx, dy);
#else
    value_t px = cosphi * dx - sinphi * dy + moffsetx;
    value_t py = sinphi * dx + cosphi * dy + moffsety;

    int lx = v2si(px * m_1dbs);
    int ly = v2si(py * m_1dbs);

    value_t x = px * m_1dbs - lx;
    value_t y = py * m_1dbs - ly;
    //printf ("%d: %.12g %.12g, lx=%ld, ly=%ld, x=%g, y=%g\n", level, px, py, lx, ly, x, y); // sanity debug

    value_t ret = 0;

    #define addon(i,j) \
          { \
            value_t bcx = 0.5 + i; \
            value_t bcy = 0.5 + j; \
            value_t dist2 = (x - bcx) * (x - bcx) + (y - bcy) * (y - bcy); \
            value_t hf = dist2hf (dist2); \
            if ( hf != 0 ) \
                ret += db.get_by_layer (lx + i, ly + j) * hf; \
          }

    addon(-1, -1);
    addon(-1, 0);
    addon(-1, 1);
    addon(0, -1);
    addon(0, 0);
    addon(0, 1);
    addon(1, -1);
    addon(1, 0);
    addon(1, 1);

    return ret * height;
#endif
  }

  inline double cov (int64_t x1, int64_t y1, int64_t x2, int64_t y2) const __PERC_INLINE
  {
    assert(0);
    return 0;
  }

private:
  void reset ()
  {
    reset_offset ();
    db.reset();
  }

  inline void reset_offset () __PERC_INLINE
  {
#ifndef UNCORRELATED
    moffsetx = 2 * m_blocksize + (SQRT2 - 1.0) * FIELD_SIZE/2 + m_blocksize * rndf.rnd1() + FIELD_SIZE / 2;
    moffsety = 2 * m_blocksize + (SQRT2 - 1.0) * FIELD_SIZE/2 + m_blocksize * rndf.rnd1() + FIELD_SIZE / 2;
    double phi = rndf.rnd1() * 2 * M_PI;
    cosphi = cos(phi);
    sinphi = sin(phi);
    //printf ("%d: moffsetx=%g, moffsety=%g\n", level, moffsetx, moffsety);
#endif
  }

  friend class Field;
};

struct pos_t
{
  typedef int c_t;
  c_t x;
  c_t y;

  inline pos_t() { }
  inline pos_t(c_t x, c_t y) : x(x), y(y) { }

  friend inline bool operator== (const pos_t& l, const pos_t& r)
  { return l.x == r.x && l.y == r.y; }
  friend inline bool operator!= (const pos_t& l, const pos_t& r)
  { return !(l == r); }

  struct hash_t
  {
    inline size_t operator() (const pos_t& pos) const
    {
      return (pos.x << 16) + pos.y;
    }
  };
};

struct coord_t
{
  struct pos_t pos;
  bool low;

  inline coord_t() { }
  inline coord_t(pos_t::c_t x, pos_t::c_t y, bool low) : pos(x, y), low(low) { }

  inline void hamster ()
  {
    if ( (pos.x ^ pos.y) & 1 )
      {
        if ( low )
            pos.y++;
        else
            pos.y--;
        low = !low;
      }
    else
      {
        if ( low )
            pos.x++;
        else
            pos.x--;
      }
  }

  inline void rhamster ()
  {
    if ( (pos.x ^ pos.y) & 1 )
      {
        if ( low )
            pos.x--;
        else
            pos.x++;
      }
    else
      {
        if ( low )
            pos.y++;
        else
            pos.y--;
        low = !low;
      }
  }

  inline void punch ()
  {
    low = !low;
  }

#if 0
  inline bool same_box (const coord_t& c) const
  {
    // vertical
    if ( pos.x == c.pos.x &&
        (pos.y == c.pos.y - 1 || pos.y == c.pos.y + 1) &&
        (low ^ c.low) )
        return true;

    // horizontal
    if ( pos.y == c.pos.y &&
        (pos.x == c.pos.x - 1 || pos.x == c.pos.x + 1) &&
        !(low ^ c.low) )
        return true;

    return false;
  }
#endif

  friend inline bool operator== (const coord_t& l, const coord_t& r)
  { return l.pos == r.pos && (l.low == r.low); }

  friend inline bool operator!= (const coord_t& l, const coord_t& r)
  { return !(l == r); }
};

class Field
{
public:
  Field (bool add_baseline = false);

  value_t value (double xc, double yc);

  inline bool valueb (int xc, int yc) __PERC_INLINE
  {
    return value (xc, yc) >= sealevel;
  }

  static inline int midpos() __PERC_INLINE { return FIELD_SIZE / 2; }

#if 0
  static void calc_cell_dist (const char *path);
  static void calc_dist (const char *path, size_t samples = 100000, size_t levels = 30);
#endif
  void printsurf (const char *path = "surf.dat", int width = 177, int height = 100, bool boolean = false);
  void testdim ();

  void reset ();
  void reset_offsets ();

  static inline bool upw (int xc, int yc) { return (xc ^ yc) & 1; }

  void set_sealevel (const value_t& sl) { sealevel = sl; }

private:
  std::vector<Layer *> layers;

  value_t sealevel;
  bool add_baseline;

#ifndef UNCORRELATED
  value_t baseline;
  value_t baselinernd;
  Database db0;
  BASERNG baserng;
#endif
};

class Walker
{
public:
  Walker (Field *field, const coord_t& start) : field (field), start(start) { }

  int prewalk (int limit = 1024);

  bool walk (int limit = 1024);
  bool walk_perim (int limit = 1024);
  bool walk_perim2 (int limit = 1024);

  typedef std::vector<coord_t> path_t;
  void print_path (const char *filen, const path_t& path);
  void print_perimpath (const char *filen, const std::vector<pos_t>& path);
  typedef std::vector<pos_t> grid_t;
  static void print_grid (const char *path, const grid_t& grid);
  void print_surf (const grid_t& grid);

  inline size_t get_hull () const { return hull; }
  inline size_t get_bonds () const { return bonds; }
  inline size_t get_ends () const { return ends; }
  inline size_t get_lines () const { return lines; }
  inline size_t get_sides () const { return sides; }
  inline size_t get_holes () const { return holes; }
  inline size_t get_perim () const { return perim; }
  inline size_t get_fjords () const { return fjords; }
  inline size_t get_corners () const { return corners; }
  inline double get_len () const { return len; }

private:
  Field *field;
  coord_t start;
  size_t hull, bonds, ends, lines, sides;

  size_t perim, holes, fjords, corners;
  double len;
};

class aggr_freq
{
  struct entry
  {
    entry() { }
    entry (int box_size, int property, int value) : box_size(box_size), property(property), value(value) { }
    int box_size;
    int property;
    int value;
    friend bool operator< (const entry& l, const entry& r)
    {
      if ( l.box_size < r.box_size )
          return true;
      if ( l.box_size > r.box_size )
          return false;
      if ( l.property < r.property )
          return true;
      if ( l.property > r.property )
          return false;
      if ( l.value < r.value )
          return true;
      return false;
    }
  };
  typedef std::map<entry,size_t> adata_t;
  adata_t adata;

public:
  void add (int box_size, int a1)
  {
    int a[] = {a1};
    dispatch(box_size, a, 1);
  }

  void add (int box_size, int a1, int a2)
  {
    int a[] = {a1, a2};
    dispatch(box_size, a, 2);
  }

  void add (int box_size, int a1, int a2, int a3)
  {
    int a[] = {a1, a2, a3};
    dispatch(box_size, a, 3);
  }

  void add (int box_size, int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8, int a9)
  {
    int a[] = {a1, a2, a3, a4, a5, a6, a7, a8, a9};
    dispatch(box_size, a, 9);
  }

  void add (int box_size, int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8, int a9, int a10, int a11, int a12, int a13, int a14, int a15, int a16, int a17, int a18)
  {
    int a[] = {a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15, a16, a17, a18};
    dispatch(box_size, a, 18);
  }

  void print(FILE *fp)
  {
    for ( adata_t::iterator it = adata.begin(); it != adata.end(); it++ )
        fprintf (fp, "%d,%d,%d,%zu\n", it->first.box_size, it->first.property, it->first.value, it->second);
  }

  void add (FILE *fp)
  {
    char line[1024];
    while ( fgets(line, 1024, fp) )
      {
        if ( !isdigit(line[0]) )
            continue;
        entry e;
        size_t count;
        if ( sscanf(line, "%d,%d,%d,%zu\n", &e.box_size, &e.property, &e.value, &count) == 4 )
            adata[e] += count;
      }
  }

private:
  void dispatch (int box_size, int* d, int count)
  {
    for ( int i = 0; i < count; i++ )
        adata[entry(box_size, i, d[i])] += 1;
  }
};

template <size_t M, bool dof = false>
class aggregator
{
public:
  inline aggregator() { clear(); }

  inline void add (int a1)
  {
    int a[1] = {a1};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (long double a1)
  {
    long double a[1] = {a1};
    for ( size_t i = 0; i < 1; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (long double a1, long double a2)
  {
    long double a[2] = {a1, a2};
    for ( size_t i = 0; i < 2; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (long double a1, long double a2, long double a3)
  {
    long double a[3] = {a1, a2, a3};
    for ( size_t i = 0; i < 3; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4)
  {
    int a[4] = {a1, a2, a3, a4};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5)
  {
    int a[5] = {a1, a2, a3, a4, a5};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6, int a7)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6, a7};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6, a7, a8};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8, int a9)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6, a7, a8, a9};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8, int a9, int a10, int a11, int a12)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8, int a9, int a10, int a11, int a12, int a13, int a14)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  inline void add (int a1, int a2, int a3, int a4, int a5, int a6, int a7, int a8, int a9,
      int a10, int a11, int a12, int a13, int a14, int a15, int a16, int a17, int a18)
  {
    int a[M] = {a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15, a16, a17, a18};
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += a[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += (long double)a[i] * a[j];
      }
    count++;
  }

  void add (const aggregator& aggr)
  {
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] += aggr.sum[i];
        for ( size_t j = i; j < M; j++ )
            cross[i][j] += aggr.cross[i][j];
      }
    count += aggr.count;
  }

  void clear ()
  {
    for ( size_t i = 0; i < M; i++ )
      {
        sum[i] = 0;
        for ( size_t j = 0; j < M; j++ )
            cross[i][j] = 0;
      }
    count = 0;
  }

  const char *sprint ()
  {
    char tmp[1024];
    str = "";
    sprintf (tmp, "%zu,%zu", M, count);
    str += tmp;
    for ( size_t i = 0; i < M; i++ )
      {
        if ( dof )
            sprintf (tmp, ",%.8Le", sum[i]);
        else
            sprintf (tmp, ",%.0Lf", sum[i]);
        str += tmp;
      }
    for ( size_t i = 0; i < M; i++ )
      {
        for ( size_t j = i; j < M; j++ )
          {
            if ( dof )
                sprintf (tmp, ",%.8Le", cross[i][j]);
            else
                sprintf (tmp, ",%.0Lf", cross[i][j]);
            str += tmp;
          }
      }
    return str.c_str();
  }

  inline bool empty () const { return count == 0; }

  long double avg(int n) const { return sum[n] / count; }
  long double var(int n) const { return cross[n][n] / count - avg(n) * avg(n); }
  size_t size () const { return count; }

private:
  long double sum[M];
  long double cross[M][M];
  size_t count;
  std::string str;
};

class freq
{
public:
  inline void add (int L) __PERC_INLINE
  {
    f[L]++;
  }

  void print (const char *basename, int A);

private:
  typedef std::map<int,int> f_t;
  f_t f;
};

